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HIGHLIGHTS 


• A numerical thermodynamic model for Stirling engine systems was developed. 

• Thermodynamic equations were coupled with the heat transfer governing equations. 

• The model was validated with experimental and numerical data. 

• The brake power and engine efficiency at different conditions were calculated. 

• Additional model results provide a deeper insight into the engine operation. 
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The reliability of modelling and simulation of energy systems strongly depends on the prediction ac¬ 
curacy of each system component. This is the case of Stirling engine-based systems, where an accurate 
modelling of the engine performance is very important to understand the overall system behaviour. In 
this sense, many Stirling engine analyses with different approaches have been already developed. 
However, there is a lack of Stirling engine models suitable for the integration into overall system sim¬ 
ulations. In this context, this paper aims to develop a rigorous Stirling engine model that could be easily 
integrated into combined heat and power schemes for the overall techno-economic analysis of these 
systems. The model developed considers a Stirling engine with adiabatic working spaces, isothermal heat 
exchangers, dead volumes, and imperfect regeneration. Additionally, it considers mechanical pumping 
losses due to friction, limited heat transfer and thermal losses on the heat exchangers. The model is 
suitable for different engine configurations (alpha beta and gamma engines). It was developed using 
Aspen Custom Modeller® (ACM®) as modelling software. The set of equations were solved using ACM® 
equation solver for steady-state operation. However, due to the dynamic behaviour of the cycle, a C++ 
code was integrated to solve iteratively a set of differential equations. This resulted in a cyclic steady- 
state model that calculates the power output and thermal requirements of the system. The predicted 
efficiency and power output were compared with the numerical model and the experimental work 
reported by the NASA Lewis Research Centre for the GPU-3 Stirling engine. This showed average absolute 
errors around ±4% for the brake power, and ±5% for the brake efficiency at different frequencies. 
However, the model also showed large errors (±15%) for these calculations at higher frequencies and low 
pressures. Additional results include the calculation of the cyclic expansion and compression work; the 
pressure drop and heat flow through the heat exchangers; the conductive, shuttle effect and regenerator 
thermal losses; the temperature and mass flow distribution along the system; and the power output and 
efficiency of the engine. These results show that the model allows an extensive study of different pa¬ 
rameters of the engine and thus it is suitable for design optimization studies. In addition, it also presents 
the capability for the integration into overall Stirling engine combined heat and power systems and 
therefore will allow the performance evaluation of the engine integrated on these systems. 
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Nomenclature 

St 

non-dimensional Stanton number 


area (m 2 ) 

SE 

Stirling engine 

A 

T 

temperature (K) 

A 0 

external wet area of the tube (m 2 ) 

Tad 

adiabatic flame temperature of the fuel (K) 

c f 

non-dimensional friction coefficient 

Twi 

temperature at the internal wall of the tubes (K) 

Cfd 

form drag coefficient 

Two 

temperature at the outer wall of the tubes (K) 

C sf 

skin friction coefficient 

Twater_in 

inlet temperature of the water (K) 

Cp 

constant pressure specific heat (J/kg K) 

U 

mean velocity (m/s) 

Cpwater 

constant pressure specific heat for inlet water (J/kg K) 

V 

mean velocity (m/s) 

D d 

displacer diameter (m) 

V 

volume (m 3 ) 

di 

internal diameter of the tube (m) 

W c 

compression work (J/cycle) 

dhy 

hydraulic diameter (m) 

W e 

expansion work (J/cycle) 

Err 

error tolerance 

W e , 

electric work output (W) 

Errorl 

absolute error calculated for T c and T e 

w n 

net work by cycle (J/cycle) 

Error2 

absolute error calculated for T k and T h 

Wploss 

energy loss due to pressure drop (J/cycle) 

Error3 

absolute error calculated for T wh and T wk 

W br 

engine net brake work (J/cycle) 

f 

friction factor coefficient 

A P 

pressure drop (N/m 2 ) 

freq 

engine frequency (1/s) 

z 

displacer stroke (m) 

Fr 

view factor 



h 

convective heat transfer coefficient (W/m 2 I<) 

Acronyms 

hr 

radiation heat transfer coefficient (W/m 2 I<) 

ACM 

Aspen custom modeller 

hwater 

water film heat transfer coefficient (W/m 2 I<) 

CHP 

Combined heat and power 

J 

annular gap cylinder displacer (m) 

LeRC 

Lewis Research Center 

k 

thermal conductivity (W/m K) 



L 

tube length (m) 

Subscripts 

m 

mass of the fluid (kg) 

c 

compression space 

M 

molecular weight (kg/mol) 

d 

displacer 

141 water 

mass flow of the inlet water (kg/s) 

k 

cooler space 

NTU 

number of transfer units 

e 

expansion space 

Nu 

non-dimensional Nusselt number 

f 

final value 

P 

pressure (Pa) 

h 

heater space 

Pbr 

engine brake power (W) 

r 

regenerator space 

Pr 

non dimensional Prandtl number 

i 

inside section 

Q 

heat transfer rate (W) 

0 

outside section 

Qhc 

heater heat transfer rate by cycle (J/cycle) 

w 

wall 

Qkc 

cooler heat transfer rate by cycle (J/cycle) 

wg 

wetted gas 

Qrc 

regenerator heat transfer rate by cycle (J/cycle) 

pist 

piston 

Qht 

total heating requirement for the engine (W) 

0 

initial value 

Qkt 

total cooling requirement for the engine (W) 



Qlossr 

heat loss due to imperfect regenerator (W) 

Greek symbols 

Qlk 

heat loss due to internal conduction (W) 

a 

surface emissivity 

Qlsh 

heat loss due to shuttle conduction (W) 

T 

adiabatic constant 

R 

gas constant (J/kg K) 

Vb 

brake efficiency 

Re 

non-dimensional Reynolds number 

a 

Stefan-Boltzmann constant (W/m 2 K 4 ) 

Rci 

conductive thermal resistance for tubes wall (K/W) 

e 

regenerator effectiveness 

Rfi 

fouling thermal resistance inside the tubes (K/W) 

P 

fluid density (kg/m 3 ) 

Rf 0 

fouling thermal resistance outside the tubes (K/W) 

0 

crank angle (rad) 

Rhi 

convective thermal resistance inside the tubes (K/W) 


viscosity (kg/m s) 


1. Introduction 

Increasing energy demands, strict environmental regulations 
and availability of different renewable resources requires the 
development of environmentally friendly energy systems. This 
needs an intensive research for the development of efficient and 
sustainable power technologies. In this scenario, fuel flexible and 
energy efficient technologies, such as the Stirling engine (SE), are 
becoming renewed solutions to address these requirements. 

Research on these engines was extensive since its invention but 
the fast growth of parallel technologies, such as internal 


combustion engines, slowed its progress. However, the develop¬ 
ment of new materials combined with engineering advances and 
the potential advantages that the Stirling engine offers in terms of 
low emissions, low noise and fuel flexibility, has renewed the in¬ 
terest on its development [1—3]. In addition, the successfully 
integration with micro and small scale combined heat and power 
(CHP) systems [3-6] increased its potential to become a strong 
technological alternative. However, Stirling-CHP technology must 
overcome practical limitations, such as the lower electrical per¬ 
formance compared with other technologies [7 . These limitations 
arise from the complexities of the system, which must considers 
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the engine mechanical arrangement, the heat transfer through the 
different components, the fluid dynamics of the working gas in the 
engine and an adequate power control system. Consequently, the 
development of CHP systems needs a complete analysis of the in¬ 
tegrated engine performance considering different geometrical and 
thermodynamic parameters. Therefore, the combination of exper¬ 
imental and simulation tests are essential for its design, evaluation 
and optimization. For this reason, the development of a model tool 
that simulates in detail the Stirling engine performance and could 
be easily integrated into combined heat and power schemes for an 
overall analysis of these systems is proposed in this paper. 

The ideal cycle analysis is the simplest SE thermodynamic 
model. However, real performances largely differ from this ideal 
calculation. Considering this, different mathematical models have 
been developed to predict actual engine performances. According 
to Dyson et al. [8], these models could be categorized into zero 
order, approximate (first order), decoupled (second order), nodal 
analysis (third order), and the recently multidimensional analysis 
(fourth order). 

Zero order methods estimate the power output considering the 
overall size of the engine. Beale [9] presented a correlation that 
estimates the power output of the engine and his approach was 
later completed by Senft 10], who derived the so called generalized 
Beale number. These methods are empirical and thus predict 
inaccurate power outputs, only suitable for engines under certain 
geometrical and operational conditions. First order methods esti¬ 
mate the power output and the work and heat flow during the 
engine cycle. These are based on Schmidt 8] analysis, who obtained 
closed-form solutions for the special case of sinusoidal volume 
variations and isothermal hot and cold spaces. Consequently the 
model accuracy is increased, but the isothermal assumption is still 
an idealized case. The modelling accuracy is improved with the 
second order or decoupled methods, which begin with a simplified 
cycle analysis for the estimation of the cyclic heat and work. Then, 
the power and heat losses, which are assumed independent, are 
evaluated and subtracted from this first approximation. Recent 
works such as Formosa 11 , Cheng and Yu [12 , Mehdizadeh and 
Stouffs [13 , Urielli [14], and Parlak et al. [15 are based on this 
decoupled approach. This approach presents a reasonable accuracy 
with a fast solution, and thus became ideal for first design stages 
and optimization studies. On the other hand, third order nodal 
techniques increased the level of complexity and the accuracy in 
the models. These methods divide the engine into a network of 
nodes, and setup differential equations based on the conservation 
of mass, momentum and energy. Therefore, the decoupled 
assumption is replaced by considering the interaction between the 
different energy losses. Some of the analyses at this level were 
developed by Finkelstein 16], Urielli 17], Schock [18 , Gedeon [19] 
and Tew et al. [20]. In addition Organ [21 and Larson [22], 
simplified the model solution using the method of characteristics 
for the case of unsteady one-dimensional flow. These methods 
present a detailed description of the engine performance, consid¬ 
ering the heat and work flows and the thermal and mechanical 
losses. However, the computational time is largely increased and 
thus it is not suitable for optimization studies. Finally, fourth order 
methods, or multidimensional analysis has been recently devel¬ 
oped thanks to the use of computer fluid dynamic (CFD) simulation 
codes. These multi-dimensional models provide detailed informa¬ 
tion regarding the flow pattern, temperature, and pressure distri¬ 
butions inside the engine. But, it requires much more computation 
efforts which are expensive and restricted to detailed engine 
design. This approach was used by Mahkamov [23 , Ibrahim [24] 
and Wilson et al. [25]. 

The modelling techniques that were previously described are 
focused on stand-alone engine analysis and the suitability for the 


integration into overall systems is not considered. Few works 
studied in detail this integration by simulation techniques 
[26-30]. These works are mainly based on the Annex 42 empir¬ 
ical model [31,32], which was developed for a particular engine 
and thus needs experimental data to adjust the results for other 
engine configurations. For this reason, this approach is closer to a 
“black box model” that does not reflect the thermodynamics of 
the engine. Therefore, the objective of this work is to develop an 
alternative Stirling engine model, which describes the thermo¬ 
dynamics, heat transfer and fluid dynamics through the engine. In 
addition, the modular approach considered smoothes the inte¬ 
gration within overall CHP schemes and thus it is suitable for 
different engine configurations. Thereby, this model reflects the 
engine thermodynamics and it is suitable for overall CHP techno- 
economic analysis. 

Aspen Custom Modeller (ACM) was chosen as modelling tool 
due to its modular flow sheeting environment, algebraic differen¬ 
tial equations solver, dynamic systems simulation and the possi¬ 
bility to link with external code sources (C++). The model 
proposed is based on the work of Urielli 14], assuming adiabatic 
working spaces and isothermal heat exchangers. The energy losses 
due dead volumes, pressure drop, limited heat transfer and 
imperfect regeneration are included. Additionally, external heat 
balances for the engine integration with the heat sources in the hot 
and cold side of the heat exchangers is considered. This model is 
applied for the simulation of the GPU-3 Stirling engine integrated 
within the external burner and a water cooling system as described 
on LeRC reports [33]. 

2. Second order model numerical formulation 

The model for the Stirling engine system consists of four main 
modules named ideal adiabatic, internal heat transfer, external heat 
transfer, and energy losses. The system scheme and the relationship 
between the modules are shown in Fig. 1. 

The first module represents the second order Stirling engine 
model developed by Urielli [14 .This model assumes ideal adiabatic 
compression and expansion spaces to estimate the main engine 
energy variables. The outputs are coupled to the internal heat 
transfer module, which through appropriate correlations evaluates 
the heat transfer and the temperature of the working fluid inside 
the heat exchangers. Fig. 1 shows the relation between both mod¬ 
ules and the loop represents the iterative calculation of the working 
fluid temperature. 

The next module includes the heat transfer between the 
external walls, considering the heat and cold sources. This is done 
through energy balances and heat transfer correlations, similarly to 
Hsu's approach [34 .There is also an iteration loop because the 
calculated external wall temperatures were approximated inputs 
for the adiabatic module. Finally, the last module includes energy 
losses to achieve results closer to the real engine performance. 

The main variables that connect the modules are described 
below. Additionally, the block diagram shown in Fig. 1 is expanded 
into the scheme in Fig. 2, to represent these variables. This reflects 
the modular approach proposed to model the system: 

- External heat transfer module. This module considers the 
adiabatic flame temperature and the inlet temperature of the 
cooling fluid on the hot and cold side respectively. Therefore, the 
heat source (Qh) and the heat sink (Qk) are used to estimate the 
wall temperatures (T wo h, T wo k). This approach is proposed to 
couple the Stirling engine within the external heat and cooling 
sources respectively. 

- Internal Heat Transfer module. The internal working gas 
temperatures (Th, Tj<) in the heater and cooler respectively are 
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Design variables 



Fig. 1. Block diagram for the Stirling engine system. 


calculated using heat transfer correlations for steady-state in¬ 
ternal forced convective flow [35 . On the other hand, the 
regenerator analysis proposes the use of cyclic flow heat transfer 
correlations, which are more suitable for the flow conditions on 
this space [36 . Therefore, with these correlations the effect of 
limited heat transfer inside the engine is introduced on the 
model. 

- Ideal Adiabatic module. The main operative variables such as 
electrical output (W e i), heat and cooling demands (Qh, Qk), are 
calculated considering the internal working fluid, temperature 
distribution, and the engine geometric characteristics following 
Urielli 14] approach. 

- Energy losses Module. The losses inside the engine are esti¬ 
mated to correct the ideal adiabatic outputs. This module con¬ 
siders the losses due to pressure drop, axial conduction, shuttle 
heat transfer, and imperfect regeneration. 


2.1. Stirling engine ideal adiabatic module 

The model is based on Urielli 14] assumptions. The engine was 
divided into five different cells or control volumes that correspond 
to the expansion, heater, cooler, regenerator and the compression 
spaces. Additionally, each cell is described by its corresponding 
state variables (volume, mass, temperature and pressure), and flow 
variables (heat, work and mass flow). 

The set of equations for the adiabatic model is obtained through 
corresponding mass and energy balances applied to each control 
volume. The complete set of algebraic differential equations is 
given by Urielli [14] and also in Appendix A. This set consists of 
twenty two variables and sixteen derivatives to be solved over a 
complete cycle. Considering this, Aspen Custom Modeller was 
coupled with an external C++ code, to solve the problem iteratively 
until cyclic steady-state conditions are reached. 



Fig. 2. Scheme of the Stirling engine model. 
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2.2. Heat transfer model 

The heat transfer model considers the thermal energy that flows 
in and out of the engine. The heat flow is transferred to the working 
fluid via the heat exchangers at the hot and cold sides. This is shown 

in Fig. 3. 

The heat transfer model is divided into: 

- The external heat transfer module, which considers the heat 
flow from the hot or cold source to the external walls of the heat 
exchangers. This model considers that the main heat transfer 
mechanism from the hot source to the heater wall is by radia¬ 
tion. This assumption was probe feasible at higher tempera¬ 
tures, as it is shown in the work of Hsu et al. 37]. 

- The internal heat transfer module which considers the heat flow 
from the heat exchanger external walls, to the working fluid 
inside the heat exchanger 

Heater 



2.2.1. Internal heat transfer module 

The ideal heat transfer assumption in the adiabatic model, 
implied an infinite heat transfer coefficient (hi), and consequently 
the temperature of the heat exchangers walls were equal to the 
working gas temperatures (Th, Tk). However, this assumption is 
very ideal and the heat transfer characteristics considering con¬ 
duction, convection and fouling in both the heater and cooler heat 
exchangers are included in this module. 


2.2. 1.1. Heater. The heat transfer phenomenon from the heater wall 
to the working gas is considered as: Conductive heat transfer 
through the outer to the inner wall surface, and convective heat 
transfer from the inner wall surface to the working gas. Addition¬ 
ally, according to Hsu [37], it is considered that the heat flux is 
constant, normal to the surface and the radiation and axial heat 
transfer mechanisms are neglected. 


Cooler 



Fig. 3. Heat transfer model. 
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The conduction heat transfer through the wall follows the 
Fourier law of conduction and the convection heat transfer rate is 
consistent with Newton's cooling Law. This is expressed in terms of 
thermal resistances. Thus, considering an additional fouling resis¬ 
tance for the internal wall of the tubes (Rfih), the heat flux is 
expressed with the following equation. 



1 


Rcih + Rhih + R 


hih 


Tih 


Twoh Th) 



Where, R C ih, Rhih, Rfih Are the conductive, convective and fouling 
thermal resistances respectively. The conductive resistance de¬ 
pends mainly on the heat exchangers thermal properties. On the 
other hand the convective resistance mainly depends on the flow 
characteristics which are directly related to the convection coeffi¬ 
cient (hih). Therefore considering that the working gas is pressur¬ 
ized at high density and moving with high velocity, this coefficient 
is estimated with the following expression derived from the Rey¬ 
nolds analogy [35] . 



Cfh ‘ Cph ‘ Ph 
2cl lh Pr h 



where Cf is the friction factor, which is calculated with Blasius 
Correlations [35]. It is important to mention that the transport 
properties depends on the working gas temperature at the heater 
(Th), which is assumed to be constant along the heat exchanger 


2.2.2.2. Cooler. The heat transfer goes from the working gas inside 
the cooler space to the cooler outer wall. This assumes the 
following mechanism: Convective heat transfer from the working 
gas to the inner wall surface, and conductive heat transfer from the 
inner to the outer wall surface. Considering that the assumptions 
are analogous to the heater module, the heat flux in the cooler 
space is estimated with a similar expression: 



1 

^cik ^hik + Rfii< 


Twik Tk) 



The gas flow is also pressurized and moving with high velocity 
inside the cooler, therefore the internal heat transfer coefficient is 
also estimated with the Reynolds analogy. 



Cfk' Cpk ’ Bi< 

2d ik Pr k 



2.2.2.3. Regenerator. The regenerator is a simplified thermal model 
that considers the geometrical characteristics, matrix packing, ma¬ 
trix porosity, and the influence of different materials to evaluate the 
non-ideal regeneration. It is assumed that a unidirectional and 
incompressible working fluid flows from the heater to the cooler at 
crank angles corresponding to 0-180°, and from the cooler to the 
heater at angles 180-360°. During the first blow (0-180), the heat is 
rejected from the hot working fluid to the regenerator material and 
during the reverse flow the heat is transmitted back to the working 
fluid. Considering these assumptions, the imperfect regeneration is 
evaluated in terms of the regenerator effectiveness (e), which is 
defined according to the enthalpy change [12]. Based on this defi¬ 
nition, the heat loss from the regenerator is determined by: 

Q-lossr = (1 - e) X Q r (5) 

The regenerator effectiveness could also be expressed consid¬ 
ering the following definition: 


NTU 
1 +NTU 



NTU is the number of transfer units, which is related with the 
Stanton Number (St), and the wetted area in the regenerator [38 . 


NTU = * (7) 

The Stanton number is calculated from Nusselt number corre¬ 
lations evaluated by Thomas and Pittman 36 . From these, the wire 
screens, and metal felts correlations are shown in equations (8) and 
(9) respectively. 

Nu = (l.OlO + 0.503Re° 720 ) x ((1 - 0.847) x (1 - e)) (8) 


Nu = (o.963 + 0.352Re 0794 ) x ((1 - 1.485) x (1 - e)) (9) 


2.2.2. External heat transfer module 

The external heat transfer module considers the heat flow from 
the hot or cold source to the external walls of the heat exchangers. 
This module permits the integration of the engine with different 
heat sources. In this case the integration with an external com¬ 
bustion chamber and a water cooling system is considered. 

2.2.2.2. Heater. This module considers the heat flow from the flame 
in the combustion chamber to the heater external wall. The source is 
supposed to reach the adiabatic flame temperature, and therefore 
thermal radiation is assumed as the main heat transfer mechanism. 
Considering these assumptions, the heat flow (Qh), is calculated from 
Boltzman's Law 37 ] . This is reordered in terms of thermal resistances, 
including an additional fouling resistance as shown in equation (10). 

Qh = \ p) (Tad — Twoh) (10) 

HrtAh + R f°h 

The term h r h is the radiation heat transfer coefficient, which is 
estimated with the following equation [35]: 

h r h = aoF R( T ad + T woh) 01) 

The fouling factor strongly depends on the fuel characteristics. 
For example in the case of solid bio fuels, which give high amount of 
ash residues, this factor Rf 0 h varies from 0 to 40 m 2 I</kW 39]. 

2.2.22. Cooler. The external heat transfer model at the cooler side, 
considers the heat flow from the outer surface of the cooler to the 
external cooling fluid. This fluid flows outside of the heat exchanger 
in a cross flow configuration and it is assumed that the heat transfer 
is mainly due to a convective mechanism. Additionally since water 
cooling is typically used on these systems, this model includes 
water as the main cooling fluid. However, other fluids could be 
easily implemented into the model. 

The heat flow (Qk) and the external wall temperature of the heat 
exchanger are therefore estimated combining Newton's cooling 
Law for convective heat transfer and the energy balance. 
Combining these equations, the following expression to estimate 
the external wall temperature on the cooler outer wall is obtained. 

T w ok = Twater Jn + Qk (11 a ^ 777 7 1 T ^) 

The external heat transfer coefficient (h 0 k), is estimated using 
Coulburn correlation for external flow in a bank of tubes 
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[35].Additionally a Fouling factor is assumed on the water side. 
Therefore the total heat transfer coefficient is expressed as: 



1 


+ RfnkA, 


‘water 


fok^ok 


(13) 


2.3. Energy losses 

The energy losses due pressure drop and thermal losses are 
evaluated on this module. These are assumed independent and 
calculated after reaching cyclic steady-state conditions. This com¬ 
plements the losses due dead volumes, limited heat transfer and 
imperfect regeneration considered in previous modules. 

2.3A. Pressure drop in the heat exchangers 

Tew et al. [20 , showed that combining the continuity and mo¬ 
mentum equation for one dimensional flow, and considering that 
the Stirling engine reach cyclic steady-state conditions the 
following expression, that evaluates the pressure drop over a finite 
length, is obtained. 

AP = 3-1 p v 2 L (14) 

a hy z 

The friction factor coefficient (f) is estimated using steady flow 
correlations for flow inside tubes [35 . 

2.3.2. Pressure drop in the regenerator 

The regenerator pressure drop has a significant effect on both 
efficiency and power output of the Stirling engines. According to 
Chen and Griffin 40 , the pressure drop in the regenerator accounts 
for 70-90% of the total pressure drop through the heat transfer 
components. In this sense updated correlations, evaluated by 
Thomas and Pitt [36] , which consider oscillating flow conditions are 
used in this module. Furthermore these correlations include a wide 
range of parameters such as different Reynolds numbers, Valensi 
numbers, flush ratio, wire mesh width, wire diameter and porosity. 
The following expression evaluates the pressure drop [26 : 

AP = C f i+u 2 (15) 

The pressure drop is proportional to the number of flow re¬ 
sistors (n), the square of the mean velocity (u 2 ) and the friction 
coefficient (Cf). Additionally the friction coefficient considers the 
form drag (Cfd) and skin friction (C s f), with the following 
expression: 

Cf = Cfd+ fe (16) 

Where, 

C sf = 68.556 C fd = 0.5274 wire screens 
C sf = 70.035 C fd = 0.9307 metal felts 


2.3.3. Pumping loss 

The pumping loss per cycle is calculated by integration of the 
product of the total pressure drop on the heat exchangers and the 
derivative of the expansion volume [38 . 



2.3.4. Energy losses due to internal conduction 

The temperature difference along the engine promotes an in¬ 
ternal conduction through the walls of the different components. 
Among these, it is important to consider the losses from heater to 
the cooler, due to the internal conduction through the walls of the 
regenerator 38]. These are estimated from: 



k r A r 


(Twih T wd< ) 


(18) 


Where, Qn 0 k r , L r , and A r denote the heat loss, regenerator matrix 
conductivity, regenerator length, and regenerator conductive area 
respectively. 


2.3.5. Energy losses due to shuttle conduction 

Shuttle conduction occurs due to the oscillation of the hot dis¬ 
placer across a temperature gradient. The displacer absorbs the 
heat during the hot end of its stroke and then released it during the 
cold end of its stroke. This loss is estimated with the following 
expression [41]. 

Z 2 k Dist D d , N 

Qlsh = °- 4 ^ (Te-T c ) (19) 

Where, Z is the displacer stroke, K P i S t the piston thermal conduc¬ 
tivity, D d the displacer diameter, and J the annular gap between the 
displacer and the cylinder 


2.4. Brake work and heat requirements 

The brake work is evaluated by subtracting the pumping 
(Wpioss) and mechanical losses (W mec h) to the net work calculated 
with the ideal adiabatic model. 


W br = j (PdVe + PdVc) - Wpioss - W mech (20) 

The mechanical losses are determined experimentally for the 
particular engine, through correlations that relates these losses to 
the engine frequency [42 . 

The total heat requirements are also calculated including the 
thermal losses previously evaluated. It is assumed that these losses 
are restored by the heat exchangers, increasing the total heating 
(Qht), and cooling requirements (Qkt): 

Q-ht — Q-h + Qlossr + Qjk + Qlsh Q r > 0 (21) 

04<t = Qk + (l -e).Qr Qr< 0 (22) 


3. Solution 

3.2. Model layout 

The different modules that compose the Stirling model require 
different operational and geometrical inputs. These are arranged in 
model blocks that were programmed with ACM language. Their 
interrelation is shown in Fig. 4. From this figure, it is seen that the 
heater cooler and regenerator blocks, which contain the heat ex¬ 
changers geometrical variables; and the Comp-Exp block, that 
contains the engine geometrical and operational variables; are 
directly connected to the main block called Stirling. This block 
contains the main routine, including the C++ linkage for the model 
solution, and it also requires the information from the external heat 
source, cooling fluid and working gas properties blocks. The table 
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presented in Appendix B, describes in greater detail the mentioned 
blocks. Fig. 5 

3.2. Numerical solution 

The set of algebraic differential equations for the engine and the 
heat transfer modules is solved using an iterative initial value nu¬ 
merical method. The set of equations is solved at each time step 
using a fourth order Runge Kutta scheme for the time discretiza¬ 
tion. This is repeated until cyclic steady-state conditions, which is 
reached when the difference between the assumed initial values 
and the values calculated at the end of the cycle are lower than a 
defined error. 

The numerical solution was implemented in a C++ procedure 
that is coupled to the ACM model. The following scheme describes 
the iterative steps for the solution. Fig. 5 

4. Validation of the model 

The GPU-3 Stirling engine [33] was chosen for the simulation 
and model validation. This is a single cylinder displacer engine, 
with rhombic and sliding rod seals. It is capable to produce 
approximately 7.5 kW with hydrogen working fluid at 6.9 MPa and 
3600 rpm rotational speed [33]. This engine was studied by NASA- 
Lewis Research Centre (LeRC) [33], and the geometrical dimensions 


and operational characteristics are well documented. In addition, 
the NASA-Lewis Research Centre developed a computational model 
[20], which was also compared with the ACM model predictions. 

4.1. Brake power comparison with experimental reports 

The engine brake power is defined as the net brake work per 
cycle times the engine frequency (freq). 

Pbr=W 5r xfreq (23) 

The ACM model calculates the brake power considering the 
mechanical losses measured experimentally by LeRC [43]. Table 1, 
presents the brake power predicted by the model at different fre¬ 
quencies, considering hydrogen as working fluid at a mean pressure 
equal to 2.76 MPa, and the heater temperature (Th) equal to 704 K. 
These results are compared with LeRC experimental and numerical 
data. In addition, the results for the simulation at three different 
mean pressures (P = 1.38 MPa; P = 2.76 MPa; P = 4.14 MPa) are 
shown in Fig. 6. 

The previous results showed the capability of the model under 
different pressure levels. These results are complemented with the 
brake power calculation considering different temperatures on the 
hot side of the engine. This is shown in Fig. 7. 

Table 1, and Figs. 6 and 7 show that both the ACM and LeRC 
models present similar degree of accuracy when comparing with 



Fig. 4. Stirling engine model ACM layout. 
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Fig. 5. Numerical solution for the Stirling engine system. 
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Table 1 

Estimated and predicted brake power for the GPU3 engine. 


Frequency 

(Hz) 

Estimated indicated 
power (kW) 

Estimated pumping 
losses (kW) 

Mechanical 
losses (kW) 

Brake power 
ACM (kW) 

Brake power 
LeRC (kW) 

Measured brake 
power (kW) 

Absolute error 
ACM (%) 

Absolute error 
LeRC (%) 

16.67 

1.60 

0.06 

0.48 

1.06 

1.14 

1.13 

6.15 

0.9 

25 

2.19 

0.14 

0.58 

1.48 

1.60 

1.49 

0.84 

7.38 

33.33 

2.90 

0.26 

0.67 

1.96 

2.13 

1.95 

0.84 

9.23 

41.67 

3.57 

0.43 

0.77 

2.37 

2.50 

2.39 

0.68 

4.61 

50 

4.24 

0.65 

0.88 

2.73 

2.75 

2.61 

4.83 

5.36 

58.33 

4.89 

0.92 

0.98 

2.93 

2.80 

2.70 

8.40 

3.70 


the experimental measures at different frequencies, pressure levels 
and heater temperatures. For the higher temperature (Th = 704 °C) 
and intermediate pressure level, the ACM model presents more 
accurate results than the LeRC model, with a maximum error 
around ±5% for middle frequencies. However, it tends to over 
predict the brake power at higher frequencies, within an error 
around ±8%. This pattern is similar at the low pressure operation, 
where the ACM model still presents better accuracy at lower fre¬ 
quencies, but slightly over predicts the brake power at higher fre¬ 
quencies. In the case of the high pressure operation, only one 
experimental point was reported [43 , and for this point both the 
ACM and LeRC models present similar accuracy. The brake power 
prediction at lower temperatures is also similar for both models. 
However, the LeRC model presents slightly better results at this 
temperature. 

In general the ACM and LeRC models are limited by the un¬ 
derestimation of the pumping losses through the engine. The 
use of correlations for oscillating flow conditions on the ACM 
model pretended to obtain better estimates of these losses, 
which is reflected on the accuracy at middle frequencies. On the 
other hand the use of several control volumes on the LeRC model 
captures better the variation of the gas properties through the 
engine, and thus the results at lower temperatures and higher 
frequencies were more accurate. In addition, the experimental 
report pointed out that the underestimation of the losses may be 
also due to oil contamination on the regenerator and methane 
production by thermal cracking of this contaminant. This may 


explain the higher differences at higher temperatures and fre¬ 
quencies. Therefore, the ACM model capability for the prediction 
of the brake power is good enough for the engine analysis. 
Furthermore, the second order approach requires less 
computing resources than the LeRC approach, thanks to the 
use of less number of control volumes, and thus makes the 
ACM model suitable for larger studies such as overall 
optimizations. 

4.2. Efficiency 

The engine brake efficiency is defined as the brake power 
divided by the heat flow to the engine. 

* - fc (m) 

Table 2 compares the brake power predicted by the ACM and 
LeRC models at different frequencies, considering hydrogen as 
working fluid at a mean pressure equal to 2.76 MPa, and the heater 
temperature equal to 704 °C. These are compared with the exper¬ 
imental results reported by the LeRC GPU3 tests [43]. In addition, 
Fig. 8 illustrates the variation of the engine efficiency at different 
pressure levels. 

The results reflect the difficulty to predict the brake efficiency 
for the system. This is due to the complexity to evaluate the real 
heat transfer into the engine. It is important to mention that the 
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experimental values for the heat into the engine were estimated 
considering a hot energy balance and a complete combustion of 
the fuel. Therefore the experimental results have an intrinsic 
error reflected on some spread values. Despite this difficulty, the 
ACM model presents acceptable predictions at the different 
pressure levels. The errors on the estimation are around ±8% for 
frequencies higher than 16.67 Hz, with a maximum of ±14% at 
this frequency. In addition, as seen on Fig. 8, these results are 
more accurate than the LeRC results which overestimate the 
predictions at the different pressure levels. The better capability 
that the ACM presents to predict the brake efficiency may be 
explained with the use of correlations suitable for oscillating flow 
conditions, for the heat transfer through the regenerator [43]. 
From these comparisons, it can be concluded that the capability 
for brake efficiency predictions that the ACM model presents is 
valid and adequate for the design and optimization studies that 
aims to assist. 

5. Results and discussion 

The previous section analysed the model capacity to predict the 
brake power and efficiency of the engine. However, a broader 
analysis of the engine could be addressed analysing additional re¬ 
sults that the model predicts. Therefore, this section presents 
additional outcomes for the simulation of the GPU-3 Stirling en¬ 
gine, considering the operational conditions shown in Table 3. 


Table 2 

Comparisons for the engine brake efficiency. 


Frequency 

(Hz) 

Brake 
efficiency 
ACM (%) 

Brake 
efficiency 
LeRC (%) 

Measured 
brake efficiency 

(%) 

Error 
ACM (%) 

Error 

LeRC (%) 

16.67 

20.8 

28.09 

24.10 

13.69 

16.55 

25 

23.24 

29.81 

23.45 

0.89 

27.12 

33.33 

23.93 

30.71 

25.82 

7.32 

18.94 

41.67 

24.49 

29.83 

25.64 

4.49 

16.34 

50 

24.98 

28.00 

25.68 

2.73 

9.03 

58.33 

24.30 

25.00 

23.06 

5.38 

0.24 


Furthermore, some of these results are compared with the LeRC 
model calculations. 

5.1. Volume variation 

Fig. 9, shows the ACM model and LeRC results for the volume 
variation in the expansion and compression spaces during a com¬ 
plete cycle. The curves are sine waves corresponding to the dy¬ 
namics of the piston and displacer connected through a rhombic 
drive crank mechanism. These variations were determined from 
the geometric relations among the parts of the mechanism [20]. 
The mechanism provides a symmetric motion between both piston 
and displacer for this configuration. 

5.2. Pressure variation 

It was assumed that the pressure is uniform along the engine 
but it has a transient variation during the cycle. These variations are 
shown in Fig. 10 and from these a slight difference between the 
ACM and LeRC results are noticed. The difference is explained 
because the LeRC model includes a coupled correction for the 
pressures losses along the heat exchangers. On the other hand, the 
ACM model considers decoupled pressure drops, and thus these are 
subtracted after the steady-state condition is reached. Despite of 
these differences, both approaches computed a similar sinusoidal 
variation during the cycle but the ACM approach requires less 
computing effort. 

Fig. 11 shows the variation of the volumes inside the expan¬ 
sion and compression spaces and the pressure during the cycle. 
These pressure-volume plots were computed with the ACM 
model and aimed to reflect the behaviour of each working space. 
The integral of each curve correspond to the expansion and 
compression work shown in Table 4. From these curves it can be 
noticed that the net work may be increased with higher pressures 
during the expansion, but lower pressures for the compression 
processes. This will shrink the compression curve and enlarge the 
expansion one. 
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Fig. 8. Engine Efficiency at Mean Pressure = 2.76 MPa; T water in = 15 °C. 


Table 3 

Ground power unit ACM simulation run. 


Working fluid Hydrogen 

Frequency, Hz 50 

Load pressure, N/m 2 5.5 x 10 6 

Average heater temperature, K 976 

Average cooler temperature, K 300 


5.3. Mass distribution 

The mass distribution during the cycle obtained with the ACM 
model, is shown in Fig. 12. This figure shows that during the first 
part of the cycle the mass in the compression space increases, while 
the mass in the expansion space starts to decrease. This is then 



followed by a period, around t = 0.008 to 0.012, where the changes 
are less pronounced and finally a phase where the mass in the 
expansion and compression spaces noticeably increase and 
decrease respectively. 

The changes described correspond to the compression, heating¬ 
cooling, and expansion processes respectively and agree with the 
cyclic volume variations. In addition, the mass distribution through 
the other heat exchangers presents slight variations and thus re¬ 
flects the validity of using the average steady-state correlation for 
their calculations. On the other hand, it is important to notice that 
the mass on the regenerator space is lower than in the other spaces. 
This means a low pressure drop but also a poorer heat transfer rate. 
Therefore, the effect of increasing the size of the regenerator on the 
engine performance should be further studied. 


„ x 10 B 



Fig. 9. Expansion-compression space volumes. 


Fig. 10. Pressure variation as function of time. 
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Fig. 11. ACM expansion and compression pressure-volume diagrams. 


Table 4 

Engine work per cycle. 


Work per cycle ACM 


Expansion work (W e , J/cycle) 395 

Compression work (Wi c , J/cycle) 170 

Lost work due to pressure drop(J/cycle) 16 

Net work (W n> J/cycle) 209 


5.4. Heat flow 

Table 5, presents the results for the heat flow through the heater, 
cooler and regenerator at the end of a single cycle obtained with the 
ACM model. The heat flow in the regenerator is close to the ex¬ 
pected ideal zero condition at the end of the cycle. 

The cyclic variation for the heat flow is additionally shown in 
Fig. 13. The first half of the cycle is characterized by the increasing 
heat requirement on the heater, and the also increasing heat 
rejected to the regenerator. On the other hand, the last half of the 
cycle shows the heat rejected in the cooler, and the positive heat 
flow from the regenerator to the working fluid. These variations 


x IQ'* 



Fig. 12. Mass distribution along the engine obtained with the ACM model. 


Table 5 

ACM model -Stirling Engine Heat Flow for cycle. 


Heat flow 

Heater heat flow (Qi< c , J/cycle) 

395 

Cooler heat flow(Qh Cj J/cycle) 

-170 

Regenerator heat flow(Qr C J/cycle) 

-0.04 


correspond with the mass distribution and volume changes in the 
engine. Furthermore, it is important to notice the key importance of 
the regenerator, since it handles large part of the heat flow in the 
system. 

The ACM model also includes the calculation of the losses due 
to: non-ideal regeneration, conductive and shuttle effect losses. The 
variation of these losses at different engine frequencies and the 
corresponding total heat requirements are shown in Table 6. 

As seen on the table above, the changes on the internal and 
shuttle conduction losses with the operational frequency are 
almost negligible. But on the other hand, the losses due to imper¬ 
fect regeneration increases largely and thus becomes the most 
important lost at higher frequencies. It is important to consider this 
when operating at high frequencies, since the engine performance 
could be seriously affected with an inefficient regeneration process. 

5.5. Temperature variation 

The temperature inputs and fouling factors assumed for the 
calculations are shown in able 7. 

Considering these inputs, the temperatures at the heat ex¬ 
changers walls were estimated. This allowed to obtain the distri¬ 
bution shown in Fig. 14, which shows considerable differences 
between the temperature at the heat exchangers walls (T W k,T W h), 
and the temperatures of the working fluid at the different engine 
spaces(Tk, Th, T c , T e ). This reflects the correction for limited heat 
transfer included through the correlations in the model. Conse¬ 
quently, the engine efficiency computed with the model will tend 
to be lower than ideal second order approaches. In addition, 
including the fouling on the heat exchangers would further reduce 
this efficiency. However, the fouling factor should be carefully 
evaluated considering the particularities of the case studied. 

6. Conclusions 

In this paper, a Stirling engine model was developed using a 
modular methodology. This model included four modules, the 
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Table 6 

Thermal losses and requirements for the GPU3 engine obtained with the ACM model. 


Frequency 

(Hz) 

Heater flow 
(Qh. kW) 

Internal conduction 
losses (Qi k , kW) 

Shuttle conduction 
losses (Qish. 1<W) 

Regenerator 
losses (Qiossr, kW) 

Total heat requirement 
(Qh,. kW) 

16.67 

2.98 

0.47 

0.62 

0.14 

4.21 

25.00 

4.48 

0.47 

0.62 

0.26 

5.83 

33.33 

5.95 

0.47 

0.62 

0.38 

7.42 

41.67 

7.40 

0.48 

0.60 

0.52 

9.00 

50.00 

8.85 

0.48 

0.61 

0.66 

10.60 

58.33 

10.28 

0.48 

0.60 

0.80 

12.16 


Table 7 

Temperature and fouling factors. 


Tad 

Twater_in 

Rfoh 

Rfok 


1410 K 
298 K 
0 K/W 
0 K/W 



Fig. 14. Temperature variation along the engine. 


phenomena on the system. This novel approach allows a smooth 
integration of the engine model with combined heat and power 
systems and thus will allow overall design-optimization studies. The 
model also included updated heat transfer and pressure drop cor¬ 
relations for the engine heat exchangers, with a special focus on the 
regenerator that considered correlations for cyclic flow conditions. 

The proposed model was validated with experimental and nu¬ 
merical data reported for the GPU-3 Stirling engine. The results of 
the validation showed that the model present a good capability to 
compute the power output and engine efficiency under different 
frequencies, pressures and hot side temperatures. Furthermore, 
additional results for the variations of the volume, pressure, tem¬ 
perature, mass and heat flow, and energy losses through the cycle 
provided a deeper insight into the engine operation. 

Therefore, despite of some limitations the model represents an 
adequate fast tool, useful for first design and subsequent optimi¬ 
zation of the Stirling engine, with a capability for the integration 
into overall systems. In this sense, future studies will focus on the 
implementation of this modular model for the design analysis of 
overall Stirling engine based CHP plants. 
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internal and external heat transfer modules, the ideal adiabatic 
module and the energy losses module. The proposed external heat 
transfer module and the modular methodology permitted a coupling 
of the Stirling engine governing equations and the heat transfer 


Appendix A 

The scheme for the ideal adiabatic model, developed by Urielli, 
and its corresponding set of equations are shown below. 




Figure A 1. Stirling engine adiabatic model scheme. 
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Mean pressure and pressure variation 
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Appendix B 

ACM- blocks description 


Block name Main inputs 


Main outputs 


Description 


Cooler—Heater di: tubes internal diameter 

de: tubes external diameter 
len: Tubes length 
num: Number of tubes in bundle 
si: Longitudinal space between tubes 
kw: wall thermal conductivity 

Regenerator kindr: Regenerator Type 

wr: Open mesh width 
dwire: Wire diameter 
lr: Regenerator length 
por: Porosity 

Com-Exp Alpha: phase angle between 

displacer and piston 
ccl: compression clearance 
cbor: cylinder bore 
crod: connecting rod length 
crnkr:crank radius 
drod: displacer rod diameter 
dd: displacer diameter 
dstr: displacer stroke 
eel: expansion clearance 
ecc: Eccentricity 
prod:piston rod diameter 
pd: piston diameter 


a: cross sectional area of one tube 
ae: external wetted area 
aw: internal wetted area 
amin: minimal free flow area 
vk: cooler void volume 

dhy: hydraulic diameter 
air:regenerator internal area 
vr: volume regenerator 

Vclc: compression space clearance volume 
Vcle: expansion space clearance volume 
Vswc: compression space swept volume 
Vswe: expansion space swept volume 


Geometrical characteristics for the 
cooler and heater heat exchangers 


Geometrical characteristic for the 
regenerator heat exchanger 


Compression and Expansion Working 
Spaces characteristics as function of 
the drive mechanism and the engine 
configuration (alpha, beta, gamma) 


(continued on next page ) 
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( continued ) 


Block name 

Main inputs 

Main outputs 

Description 

Ext-Heat 

tad : Adiabatic flame temperature during 
fuel stoichiometric combustion 
rfeh: External fouling factor on the heater 
eair: Excess air on the combustion 

tad: Adiabatic temperature for 
the combustion 

Main characteristics for the external 
heat transfer. In this case the 
combustion process is considered 

Cooling-Fluid 

twik: Cooling fluid inlet temperature 

mwik: Cooling fluid mass flow 

fek: external fouling factor in the cooler 

twok: Cooling fluid outlet temperature 

Main Characteristic for the cooling 
external fluid. In this case the 
water is assumed 

Working GAS 

typeg: working fluid used (Air. 

Helium, hydrogen,Nitrogen) 

Rgas: universal gas constant 

Cv:constant volume specific heat 

Cp: constant pressure specific heat 

Gama: adiabatic constant 

Uo:Sutherland parameter 

To:Sutherland parameter 

So:Sutherland parameter 

Physical properties for the working gas 

Stirling 

Pmean: mean pressure 

Freq: frequency 

Additionally all the outputs 
from previous modules 

Th, Tk, Tr, Tc, Te 
dph,dpk,dpr 

Qh.Qk.Qr 

Wc,We,Wtot 

Eff: Engine efficiency 

Main module calculate the different 
requirements and operational parameters 
for the Stirling engine 
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